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Introduction 

Tbe purpose of this research effort was to begin the study of the application of /ip- 
version finite elements to the numerical solution of optimal control problems. Under NAG- 
939, the hybrid M ACS YMA /FORTRAN code GENCODE was developed which utilized 
Zi-version finite elements to successfully approximate solutions to a wide class of optimal 
control problems. In that code the means for improvement of the solution was the refine- 
ment of t h** t ime -discretization mesh. With the extension to Zip-version finite elements, 
the degrees of freedom include both nodal values and extra interior values associated with 
the unknown states, co-states, and controls, the number of which depends on the order of 
the shape functions in each element. For details, see [1]. 

One possible drawback is the increased computational effort within each element re- 
quired in implementing Zip-version finite elements. We are trying to determine whether this 
computational effort is sufficiently offset by the reduction in the number of time elements 
used and improved Newton-Raphson convergence so as to be useful in solving optimal 
control problems in real time. Because certain of the element interior unknowns can be 
^liminarpH at the element level by solving a small set of nonlinear algebraic equations in 
which the nodal values are taken as given, the scheme may turn out to be especially pow- 
erful in a parallel computing environment. A different processor could be assigned to each 
element- The number of processors, strictly speaking, is not required to be any larger than 
the number of sub-regions which are free of discontinuities of any kind. 


Summary of Completed Work 

In order to acquaint Mr. Warner better with the workings of GENCODE and with 
finite dement methods in general, much of the first part of the year was spent in the study 
and modification of GENCODE. The code as described in [1] handled only two stages, 
defined here as a time intervals with distinct differential equations. Mesh refinement was 
also prescribed beforehand such that the code always began with two elements per phase, 
doubling upon successful convergence until the final desired number discretization was 

reached. 

The modifications were conducted in parallel with a similar effort at NASA Lang- 
ley Research Center to implement the Zi-version finite elements solution procedure into a 
MATLAB optimal control problem solver, VTOTS. Since the code was being expanded to 
include an arbitrary number of stages, the governing equations were rearranged to better 
t alfP advantage of the sparse matrix solvers included in the Harwell subroutine library. 
Instead of having separate initial conditions, ter minal conditions, and jump conditions at 
stage breaks, the code was modified such that all boundary conditions, include continuity 
between stages, were included in a single set of boundary conditions. 

Tbe other major advancement of GENCODE was the freedom to specify arbitrarily 
the desired time discretization in each phase. Once the solution to that discretization 
converges successfully, that solution can be interpolated in the code to provide initial 
guesses for any subsequently desired discretization. This flexibility proved most helpful in 
solving difficult problems. 
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Once completed, the new GENCODE was tested against VTOTS, with refinements 
and error correction occurring in both codes as errors were discovered. Since then, VTOTS 
has been updated to include state constraints with a more automated switching structure, 
an improvement that we plan to incorporate in the new GENCODE. 

After the GENCODE testing was completed, we formulated the structure of using 
p- version finite elements to solve two-point boundary value problems. A summary of that 
development is included below. At this stage, a time-marching algorithm has been devel- 
oped tested which utilizes p-version elements. Results were obtained for a number of 
problems, looking at accuracy and computational effort over various length tune elements, 
numbers of Gauss quadrature points, and orders of shape functions. The improvements 
in accuracy over h - version elements are dramatic, as expected, and the overall results are 
encouraging enough that we plan to move forward into the development of a two-point 
boundary problem solver. 


Higher-Order Elements 

In moving toward eventually applying p-version finite elements to the solution of 
optimal control problems, we first developed how to use them to the solve nonlinear two- 
point boundary value problems. In order to test their performance we started solving 
problems involving given initial states and only one time element (i.e., initial-value or 
time-marching problems). 


Formulation for Optimal Control Problems 

The differential equations of interest are assumed to be of the form 

x = f{x) xeR n ‘ (1) 

where / is an autonomous function of the state vector x. The boundary conditions can be 
specified at the initial time, to, the final time t/, or some combination of both. For this 
we demote 

x(to) = X 0 ^2) 

z(tf) = x f 

For solving optimal control problems, this problem formulation corresponds to the 
Euler- Lagrange equations when the control can be solved explicitly in terms of the states 
and costates from the optimality condition. Then the vector x would include both the 
states and the costates, and the boundary conditions would include the derived costate 
boundary conditions. 

The time interval is broken up into N not necessarily equal length time elements, A t,, 
such that the time at each element boundary t, is calculated as 


<1 = to 

t, =t»-i + A ti 
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i = 2, . . . , N + 1 


( 3 ) 


and we define the states at these nodes to be 


Xi= x{ti) i = 1, ...,N+ 1 


( 4 ) 


Ultimately, only the to and t f nodal values of the states are of interest as the values at the 
internal nodes will be known functions of the 'values of the states within each tune interval. 
The time within the t th element, U is expressed as U = U - 1 4- tAU where 0 < r < 1 and 


T = 


AU 


( 5 ) 


so that dt = Atidr. 

The state equations are enforced through use of Lagrange multipliers <5A inside the 
time interval and the boundary conditions are enforced as natural boundary conditions at 
the endpoints: 

[ tf 6X T [x - f{x)] dt 4- 6X t (i - *)£ = 0 (6) 

Jt o 

Substituting for the normalized time r yields 

f>t, [' [SXfii - 6XTf(z.)] dr + S X T (i-x)|;'=0 (7) 

»=1 


where the subscript denotes the element number or node for each variable. Integrating by 
parts results in 



dr 4- -,6XJj +1 xn+i = 0 


( 8 ) 


Now from Ref. 2, we define C° shape functions for the LaGrange multipliers in each element 
in terms of nodal values and internal values 

SXi = 6 Ai(l — r) + <5A i+ ir 4- 6 Ai,i (1 - t)tV& 4 f <5 A 

= 6X[ =Ati6Xi = -6Xi + 6X i+ i 4 - 6X iA (l - 2t)V6 4 - • • • + 6~X 
dr 

Here m is the order of the shape function, with m = 0 corresponding to /i-version finite 
elements. The polynomial function a m (r) is determined through a set of relationships 
developed in Ref. 2, and a' m {r) is the derivative of a m (r) with respect to r. Similarly, the 
shape functions for the values of the states internal to each element are 

Xi =Xi t i 4- £t, 2\/6 (t — 5 ) 4- • • • 4- x*, m +i/?m+i (t) 
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Here again the polynomial function j0 m (r) is a recursive relationship developed in [2]. One 
fewer 6X term is used than x terms to ultimately ensure an equal number of equations and 
unknowns. 

Substituting Eqs. (9) and (10) into Eq. (8) and taking m = 1 for simplicity results in 


^ ] J ^ j ^— SXi 4- 6Ai + i 4- tfAq^l — 2r)\/6j [xi,i + x^y/ft ( r 2 )] 

' 0 . T -l) (11) 

4- At* |$Ai(l - r) 4- 5A i+ ir 4- 6X it i(l - r)rV6J / [*<,1 + x ii2 y/6 (r - £)J j dr 
4- £AjXi — 5 A^ + i£tv+i = 0 

Integrating all the terms that do not depend on /(r) and simplifying results in 

^ J rj"* * 

iZ Ati J 0 [^A i (l-r) + 6Ai + iT + 5A i ,i(l-T)r^] / [* if i 4- x it2 V6 (t - i)J dr 

+ Y. (-6Afxi,i 4- ^A^iXi.i - 6Xj\X it2 j + SXjxi - 6\n+iXn+i = 0 

X 

Finally, grouping terms multiplying each of the LaGrange multipliers yields 

-zi,i + *i + At i [ (1 - t)/(xi) dr 
Jo J 

in, 1 — xn+ 1 + At at / rf(xs) dr 
Jo J 

-V r r l f 1 

-*i.i + *<- 1.1 + Ati J (l-T)/(xi) dr 4- Ati -1 J ^ rf(xi-x)dT 
- -i, |2 + At, J %/6(l - t)t/(x,) drj = 0 


( 12 ) 


sXT 


^N+l 


(13) 


If n is the dimension of the state vector, then Eq. (13) results in (2 N 4- l)n equations 
for (2.V + 2)n unknowns: n unknowns at each end point and 2n unknowns within each 
element. If the order of the shape function m is greater than one, then the last* block of 
Eq. (13) becomes nm equations, and x» would then depend on x itj for j = 1, . . - , m 4- 1. 
This brings the total to (N+1 + Nm)n equations and ( N(m + 1) + 2) n unknowns. Thus 
for a uni que solution to the problem, a combination of n initial and final conditions on the 
states mus t be specified. Note that when solving these equations by a Newton- Rap hson 
approach that the Jacobian matrix is block diagonal. 
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Note that if the order of the shape function m is zero, then Xi = Xi,i the /i-version 
equations result 


SAf [-xi + xi + Ati/(xi)/2] 

+6At + i |iw - ijv+i + Atyv/(^)/2] 

N 

+ [-Xi + Xi_i 4- Atj/(xi)/2 + Ati_i/(xi_i)/2] = 0 

i=2 


Time-Marching Algorithm 

To help determine the feasibility of using p-version finite elements to solve optimal 
control problems, Eqs. (13) were first incorporated into a time-marching algorithm. In this 
case, the required number of specified boundary conditions are given as initial conditions, 
i\ = x(to). Eqs. (13) then reduce to the case of only one element N = i = 1, and we 
can drop the element subscripts. Thus, the r emainin g subscripts refer only to the element 
order. Rewriting Eqs. (13) in this way one obtains 

x i = xi - At J (1 - t)/ jxi 4- x 2 V6 (t - |)j dr 
0 = -X 2 + At\/6 J (i - r )t f |xi 4- x 2 VQ (t - 5 )] dr (15) 

x 2 = xi + At J rf jxi + x 2 V6 (t - 5 )] dr 

again a. < gniming first-order shape functions. For each additional order, there is Pne more 
equation and one more unknown. 

For shape functions of arbitrary order, the first two of Eqs. (15) become m + 1 equa- 
tions, they are solved for the internal values, x, for i — 1, 2, . . . , m + 1. These values 
are then used to calculate the nodal values, X 2 from the third of Eqs. (15). These final 
values become the initial values for the next element, and the process repeats over all the 
elements. Notice that in the solution of this problem with a Newton-Raphson algorithm, 
the Jacobian matrix is, in general, fully populated. 

The process of solving these equations has been implemented in a FORTRAN sub- 
routine ii< dug a standard full-step Newton-Rapheson algorithm. The routine utilizes the 
Harwell family of subroutines for solving fully-populated linear systems. A user interface 
was developed using the symbolic manipuation software MACSYMA which generates the 
appropriate expressions for error evaluation and the Jacobian matrix. Unlike for /i- version 
finite deme nts, the quadratures in Eq. (15) in general cannot be done be done by inspection 
and similarly are too complicated for software such as MACSYMA. To perform the inte- 
grations, Gaussian quadrature with a variable number of Gauss points was implemented 
here. 
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Example Problem Formulations 

The code was run on a variety of problems for varying orders of shape functions, 
numbers of elements, and numbers of Gauss points. The goal was to determine whether 
such a code could run in a reasonable length of time for complicated problems, thereby 
determining the feasiblility of applying p- version elements to two-point boundary value 
problems. That being established, the next step was to determine the relationship be- 
tween the number of Gauss points and the accuracy of the solution, in hopes of finding a 
for mula such as in the case of /^version elements. The results of this experimenting will 
be described for a two-state linear system, two one-state systems with different types of 
nonlinear dy nami cs, and an 8-state missile system. 


Linear System 

The first problem that was examined was an idealized spring-mass system, with a 
spring constant, k and mass, m. 

mx = —kx 


with 


z(0) = 1.0 

x(0) = 1.0 


( 16 ) 


The analytical solution to this with k = m = 1 is 

x(t) = cos(t) x(l) = 0.543087528 

x(t) = — sin(t) x(l) = —0.839676091 


In Fig. 1, the log of the percentage error in the time-marching algorithm at the final 
time is plotted versus the number of Gauss points for various orders of shape functions. 
The numb er of elements in all cases was 5. The percentage error, E, was calculated in 
term*; of the exact answer from Eq. (17) as 




Vi( 1) 


( 18 ) 


where N is the number of elements, and n is the dimension of the state vector y — {x,x}. 
In subsequent examples, the error is calculated similarly with y again being the appropriate 
state vector. This error criterion gives a good idea of how errors are being propagating 
through the algorithm. 

Notice that in this case the error is minimized when the number of gauss points is one 
more than the order of the shape function in each element. This will hold true for all linear 
systems because the polynomial in r to be integrated is of no higher order than 2m + 1, 
where m is the order of the shape function. As shown in [3], using n Gauss quadrature 
points and weights results in exact integration of polynomials up to order 2n — 1. Thus 
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only m-r 1 Gauss points are needed to exactly evaluate the integrals when the order of the 
shape functions is tn. Unfortunately this result does not hold for nonlinear systems. 

Changing the number of elements adjusted the error appropriately in each case, but 
did not change the relative characteristics between the orders of the shape functions. Also, 
the number of Newton steps required was only one for all elements sizes, orders of shape 
function, and numbers of Gauss points. Thus for a linear system, the tradeoffs between 
accuracy and computational effort are very clear. 


Nonlinear System with Quadratic Nonlinearities 

The next problem complicates matters some by replacing the linear term with a 
quadratic in a first order differential equation 


x = ax 2 

x(0) = 1.0 


(19) 


where a is a constant. The analytic solution to this problem with a — 0.5 is 

*< t)= (oi^ir i(1)=20 (20) 

In Fig. 2, the log of the percentage error at the final time from Eq. (18) is again plotted 
versus the n umb er of gauss points for various orders of shape fiinctions. The results plotted 
were obtained with 5 elements, but the trends are the same for any number. In this case 
the error is again minimized for the number of Gauss points equals the order of the shape 
function plus one. What is different here is that the accuracy actually goes down as more 
Gauss po ints are added, enforcing the point described in [3] that more Gauss points does 
not always mean higher accuracy. 

The numb er of Newton steps per element decreased on average as the number of 
elements increased, but only from 6 to 5. And there was no significant decrease in the 
number of Newton steps required as the order of the element increased. Thus the high 
accuracy solutions still require considerably more computational effort. 


Nonlinear System with Infinite Order Nonlinearities 

The next problem examined replaced the quadratic nonlinearity with an exponential 
function, an infini te-order nonlinearity. 


x = e 


X 


x(0) = -1 


The analytic solution to this problem is 

x(t ) = - ln(-t + e) - x(l) = - ln(e - 1) 


( 21 ) 


( 22 ) 
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In Fig. 3, the log of the percentage error at the final time from Eq. (18) is again 
plotted versus the number of Gauss points for various orders of shape functions. The 
figure correspond to 5 elements used, but the trends are the same for any number. Here 
the line begins to blur about the optimal number of Gauss points. In all cases, the error is 
very close to the minimum with one more Gauss point than the order of the shape function, 
but depending on the order and the number of elements, one additional Gauss point can 
still improve accuracy significantly. 

It is encourag ing , however, that the number of Newton steps here did decrease with 
increased order, though only from 6 to 5 on average for each element. In this problem, 
extra iterations are rather insignificant, as the problem has only one state. Also note 
that in none of the above examples were there any cases in which the algorithm failed to 

converge. 

Nonlinear Missile Model 

To give the time-mar ching algorithm a harder test, we next tested it on an 8-state 
missile system. While the equations themselves are unimportant, it should be noted that all 
the rumlinrar system dynamics were retained, including table-lookups on highly nonlinear 
lift and drag coefficients. The missile model was integrated from launch for 5 s, when the 
engines would be throttled back. The two controls are angle of attack and roll angle, both 
of which were set nominally to 1° throughout the tune interval. 

To provide a comparison, the model was integrated using the Runge-Kutta method 
and 1000 time steps. Fig. 4 plots the error in the finite-element algorithm as calculated 
in Eq. (18) using the shooting results as the exact answer. With any fewer than 5 time 
elements, the algorithm had difficulty converging. A restricted-step Newton method [4] 
was then implemented, which alleviated some of the convergence difficulties while not 
appreciably affecting accuracy. 

Unfortunately, none of the trends we encountered before with regard to optimal num- 
ber of Gauss points held true in this case. It is clear that for 41 real life problems, the opti- 
mal number of Gauss points may have to be determined on a case-by-case basis. Therefore, 
any future algorithms will include capability for adjusting this number. Also, increasing 
the order once again brought down the number of Newton steps required for convergence. 

Future work 

The improvement, in accuracy achievable by increasing the order of the shape functions 
in all four problems presented here is exceptional. The question still remains, however, as 
to how the extra computational effort will trade off with this increased accuracy when the 
methodology is extended to two-point boundary value problems. Another open question 
is how to find a simple means for obtaining the optimal number of Gauss points for more 
complicated problems. 

Given the encour aging results we have seen using /^version finite elements so far, we 
plan to extend GENCODE to include the hp - version finite elements and begin testing it on 
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optimal control problems, including those with state constraints. There are also a number 
of theoretical aspects of these finite elements yet to be explored including the effect on 
element stability and convergence properties. 
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Fig. 1: Time-Marching Error at Final Time for Spring Mass System 



Fig. 2: Time-Marching Error at Final Time for x — x 2 
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Fig. 3: Time-Marching Error at Final Time for x = e x 
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Fig. 4: Time-Marching Error at Final Time for Missile system 
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